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Abstract 

The origin of diversification and coexistence of genes and species have been 
traditionally studied in isolated biological levels. Ecological and evolutionary 
views have focused on the mechanisms that enable or constrain species coexis- 
tence, genetic variation and the genetics of speciation, but a unified theory linking 
those approaches is still missing. Here we introduce evolutionary graphs in the 
context of neutral theories of molecular evolution and biodiversity to provide a 
framework that simultaneously addresses speciation rate and joint genetic and 
species diversities. To illuminate this question we also study two models of evo- 
lution on graphs with fitness differences, which provide insights on how genetic 
and ecological dynamics drive the speed of diversification. Neutral evolution gen- 
erates the highest speed of speciation, species richness (i.e. five times and twice 
as many species as compared to genetic and ecological graphs, respectively) and 
genetic-species diversity (i.e., twice as many as genetic and ecological graphs, 
respectively). Thus the speed of speciation, the genetic-species diversity and 
coexistence can differ dramatically depending on whether genetic factors versus 
ecological factors drive the evolution of the system. By linking molecular, sexual 
and trophic behavior at ecological and evolutionary scales, interacting graphs 
can illuminate the origin and evolution of diversity and organismal coexistence. 



20 1 Introduction 

21 One outstanding challenge in ecology and evolution is the development of an accu- 

22 rate and complete understanding of diversity across biological levels and spatial scales 

23 [21,43,59]. The neutral and nearly-neutral theories of molecular evolution [36,39,58] 

24 (hereafter NTME) were proposed to explain observations of high evolutionary rates 

25 and the maintenance of large amounts of molecular diversity within populations [37,57] . 

26 Similarly, the neutral theory of biodiversity (hereafter NTB) promises to contribute to 

27 our understanding of how species diversity is maintained in ecological systems [2,9,30]. 

28 Both theories share the same framework [28,62]. Furthermore, both theories offer a 

29 baseline from which to extend theories of evolution [4, 47] and to test the role of fre- 

30 quency and density dependent selection on the generation, evolution and maintenance 

31 of diversity at genetic and ecological levels, respectively [22,35,43,57]. 

32 Neutral models at the molecular level have considered mutation rate (/i), random 

33 fluctuations of alleles (i.e., genetic drift), and molecular constraints on producing fer- 

34 tile offspring (i.e., the genetic similarity value q^^ between any pair of individuals i 

35 and j must be higher than g™-*") as mechanisms of speciation in populations with J 

36 individuals [11,18,26]. The neutral theory of biodiversity introduces the implicit speci- 

37 ation rate at the individual level {v) where species fluctuate randomly (i.e., ecological 

38 drift) and all individuals are equivalent (i.e., neutral competitive interactions) [30,54]. 

39 Speciation is crucial to the neutral biodiversity theory (without it diversity cannot be 

40 maintained), yet the speciation parameter is simply assumed and has no basis in bio- 



41 logical processes. To integrate genetic and ecological neutral theories, we need to link 

42 the speciation rate (z/) with explicit mechanisms of speciation from neutral molecular 

43 theories [18]. 

44 Despite the striking parallels between neutral theories in population genetics and 

45 community ecology, the speed of speciation and diversity patterns at genetic and com- 

46 munity levels have rarely been studied simultaneously [3,35,41,61,63]. This raises 

47 important questions. For example, let us consider a population with J reproductive 

48 compatible individuals. This defines a completely connected graph of size J x J. Given 

49 this initial graph in a population, does neutral evolution at molecular and ecological 

50 levels speed up speciation and increase genetic-species diversity? If frequency and 

51 density dependence effects at genetic and ecological levels are important, how can we 

52 discern the speed of speciation and genetic-species diversity under neutral or natural 
63 selection scenarios? Thus, do genetic or ecological level drive the speed of speciation, 

54 genetic-species diversity and coexistence? [8,15,16,32]. 

55 In order to answer those questions we need a framework that allows us to study the 

56 molecular and ecological levels simultaneously. This framework represents an ambitious 

57 research programme - much more than can be accomplished in a single paper. Our 

58 goals here are more limited. First, we introduce evolutionary graphs [46] in the con- 

59 text of neutral theories of molecular evolution [26,36,39] and biodiversity [29,30] which 

60 suggest a promising new way to provide a general account of how neutral, positive and 

61 negative density and frequency dependent selection affect the speed of diversification 

62 and genetic-species diversity. Second, we introduce genetic and ecological graphs where 



63 the genotype-phenotype of each individual are represented as one to one or are decou- 

64 pled by the specific behavior and phenotypic plasticity of each individual, respectively. 

65 Note that in addition to the graph of reproductive individuals, we need to specify a 

66 new graph that captures the effect of the phenotypic plasticity in the system. 

67 Evolutionary neutral graphs in the context of two mutualistically interacting pop- 

68 ulations are defined as follows. Consider two randomly mating populations of size 

69 Jr and Jp where each individual has an infinitely large genome sequence subject to 

70 random neutral mutations. The initial genetic similarity values between each pair of 

71 individuals {q^ and qp) in the matrices Qr = [g]^] and Qp = [qp] are equal to 1 and 

72 mutation rates fip and fip are equal among all individuals (Fig. la). At each time 

73 step, an individual of each population is chosen for death (Fig lb). Two individuals 

74 are chosen for reproduction. Individuals have the same probability 1/Jr {1/ Jp) to be 

75 chosen for reproduction or for death (Fig Ic). The offspring of these two individuals 

76 replaced the dead individual. The newborns in R and P can be the consequence of 

77 sexual reproduction without a mutualistic interaction (i.e., facultative mutualism given 

78 hy uo > 0) or a consequence of a mutualistic interaction with individual effectiveness m 

79 between the first two chosen individuals for reproduction in community R and P (i.e., 

80 given by 1 - cj) (Fig. Id). 

81 All these elements allow us to develop models of evolution on genetic and ecological 

82 graphs with the following additions to the neutral model (Fig. 2): (1) fitness differences 

83 within each species according to the number of genetically related mating partners (i.e., 

84 genetic graphs), or to the number of trophic links with individuals in the second com- 



85 munity (i.e., ecological graphs); (2) density dependence across species, thus rare species 

86 have higher probabilities of reproduction, and (3) contrary to genetic graphs, where 

87 the offspring can inherit the high connectance from its parents increasing its reproduc- 

88 five probability, all offspring in the ecological graph start with the same reproductive 

89 probability. Let us consider first the genotype-phenotype map as one to one. In this 

90 "genotype-fitness speciation model" (hereafter GF) reproductive probabilities are set 

91 according to the total number of genetically related individuals each individual i can 

92 mate with, so that we take into account explicit fitness differences within each species 

93 (Fig. 2b). The genetic level, assuming that all the traits involved in sexual selection 

94 are under genetic control, determines the evolution of the system based on the genetic 

95 similarity among individuals. 

96 There is empirical evidence for the effect of ecological interactions mediated by 

97 phenotypic traits on the evolution of diversity [10,34,42,53,56], but they have so far 

98 been missing in neutral theories. Let us consider that the phenotype is not simply the 

99 product of the genotype, but that it is influenced by the interactions with individuals of 

100 a second community (i.e., second trophic level). In the "phenotype-fitness speciation 

101 model" (hereafter PF) we still have the genetic similarity constraint on having fertile 

102 offspring, but the role of ecological interactions is represented as a varying reproductive 

103 probability for each individual according to its specific behavior, development and 

104 phenotypic plasticity [33]. Those phenotypic traits, not linked to the total number 

105 of genetically related matings with individuals of the same species, are given by the 

106 evolution of the number of trophic links with individuals of the second community. 



107 Thus, the reproductive probabihty of each individual increases with the number of 

108 trophic hnks with individuals of the second community, but it is independent of the 

109 number of potential genetically related matings with individuals of the same species 
no (Fig. 2c). 

111 We show that the neutral scenario, which is represented by a diverse genetic pool 

112 of parents in the context of decoupled evolving mating and trophic graphs, triggers 

113 the highest speed of speciation and highest levels of genetic-species diversity and co- 

114 existence. We also show that ecological graphs, whose reproduction is determined by 

115 specific behavior or phenotypic plasticity and not by the total number of genetically 

116 related matings, allow higher speciation rate and coexistence than mating graphs. Link- 

117 ing neutral theories at the molecular and ecological levels using evolving graphs promise 

118 to contribute to our understanding of contemporary diversity at multiple levels [12]. 

119 As we will demonstrate, it represents a powerful starting point to: 1) understand the 

120 speed of speciation and the relationship between genetic and species diversity by using 

121 genetic and ecological graphs [18,26,46], and 2) understand the consequences of niche 

122 and neutral dynamics as a continuum that is based on ecological interactions among 

123 individuals [22,52,64] 

124 2 Results 

125 First, not surprisingly, mating and trophic number of links at the individual level are 

126 not correlated during the evolution of the system under the neutral and the phenotype 

127 fitness scenarios (Fig. 3a and Fig. 3b, respectively). The distribution of Spearman's 



12B rank coefficient values is close to a uniform distribution in both cases (Fig. 3a and 

129 3b represent the distribution for community Jr). The distribution of the Spearman's 

130 values for the genotype fitness model is highly skewed with approximately 80% of p- 

131 values < 0.01, suggesting that mating and trophic degree are in most cases correlated. 

132 Fig. 3d and 3e represent the evolution of individual mating and trophic degree in the 

133 genotype fitness model as a function of the individual rank (i.e., from the most (left) 

134 to the least (right) connected individual). 

135 The neutral unified model generates on average twice and five as many speciation 

136 events (i.e., 188 ± 10) as the phenotype (i.e., 80 ± 1) and the genotype fitness models 

137 (i.e., 34 ± 3) respectively (Fig. 4a, results for the community Jp not shown but with 

138 the same parameter values they are qualitatively the same). Similarly, waiting time to 

139 speciation or the number of generations to the ffist speciation event is on average twice 

140 and five times as small in the neutral case (170 ± 3) as in the phenotype (440 ± 12) 

141 and the genotype fitness scenarios, respectively (924 ± 25) (Fig. 4a, see appendix for 

142 a detailed description of the sampling of the transients and the steady state). At 

143 stationary state (approx. 1000 generations, see Fig. 5) speciation events scale with the 

144 number of generations for all the three models (r^ = 0.99) with the scaling exponent 7 

145 = 0.97 (neutral), 1.03 (phenotype fitness) and 1.31 (genotype fitness), red lines in Fig. 

146 4a. 

147 Note that we have used the same three input parameters in the three models ex- 

148 plored. Mutation rate, with /i^ = /xp = /x, the minimum genetic similarity value g™" = 

149 g™" = g™*" and the individual mutualistic effectiveness niij = rriji = m = 1 assuming 
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150 a fully symmetric case for all the individual interactions in the context of obligate mu- 

151 tualism (i.e., a; = 0) (see Methods). Thus, the speed of speciation rate is driven by the 

152 specific reproductive transition probabilities at individual level. This result remains 

153 similar after relaxing the assumptions of effectiveness and facultative or obligate mu- 

154 tualism. Does the distribution of the number of generations to speciation differ among 

155 the models? All the nontransformed distributions were highly skewed (skewness index 

156 > 2), and significantly different from a normal distribution (Fig. 4b Lilliefors's test, 

157 all P < 0.001 with means of 47, 258, and 115 for the neutral, the genotype and the 

158 phenotype scenario, respectively). The distribution of the number of generations to 

159 speciation differ significantly between the neutral and the genotype/phenotype models 

160 {Kolmogorov — Smirnov test, P < 0.0001), but the genotype and the phenotype fitness 

161 scenarios do not differ significantly {Kolmogorov — Smirnov test, P > 0.1). 

162 The neutral model generates on average twice as many genetic and species diversity 

163 as the phenotype and the genotype fitness scenarios (Fig. 5a using eq. (2) in Methods, 

164 and 5b, using species diversity S^ as ^^^^^ — j , where pi is the relative abundance of 

Z^i Pi 

165 species i). As in the speed of speciation, the neutral scenario predicts twice and 

166 five as many number of coexisting species as the phenotype and the genotype fitness 

167 model, respectively (Fig. 5c). Genetic diversity (Fig. 5a), species diversity (Fig. 5b) 

168 and species richness (Fig. 5c) distributions for all the models differ from a normal 

169 distribution {Lilliefors' test, P < 0.001) despite their strong differences in skewness. 

170 The neutral case predicts highly symmetric distributions, all skewness indices between 

171 -0.15 and 0.08, while the phenotype and the genotype model predict skewness indices 
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172 between 0.87 and 1.44 and > 2, respectively. Genetic-species diversity and species 

173 richness distributions differ significantly among all the models (Kolmogorov — Smirnov 

174 test, P < 0.0001). 

175 In summary, the diverse genetic pool underlying our unified neutral scenario in the 

176 context of the uncorrelated mating and trophic graphs triggers the highest speed of 

177 speciation with consequences to the genetic-species diversity, coexistence and species 

178 richness. Note, however, that the species diversity values with the explicit mechanisms 

179 of speciation are lower than the values from the biodiversity number 6b in the neu- 

180 tral theory of biodiversity. These results remain qualitatively similar for the range of 

181 parameter combinations explored (see appendix). 

182 3 Summary and Discussion 

183 The present study is an attempt to unify the speed of speciation with the evolution 

184 of diversity at genetic and ecological levels. We create a bridge between the neutral 

185 theory of molecular evolution [36, 39] and the neutral theory of biodiversity [30] using 

186 mating and ecological graphs in the context of explicit mechanisms of speciation [17, 

187 26]. The unified neutral model predicts the highest speed of speciation, number of 

188 coexisting species (i.e., five and twice as many as genetic and ecological networks, 

189 respectively), and genetic-species diversity (i.e., twice as many as genetic and ecological 

190 networks), but diversity values are lower than the neutral biodiversity theory with 

191 implicit speciation. This result is not surprising. Genetic variation maintained in non 

192 random mating is to same extent cryptic since the heterozygote diversity is less than 
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193 from a random mating population. However we show how the speed of speciation and 

194 genetic-species diversity are closely controlled by the dominant graph (i.e., genetic or 

195 ecological) at each level during the evolution of the system. 

196 Note that we have explored only a few scenarios (see appendix). Despite that the 

197 effect of the genetic regulatory [12,44,55] and ecological interactions [10,34,42,56] on 

198 the evolution of diversity is widely recognized, they have so far been missing in neutral 

199 theories. Here we show that the decoupling of phenotypic (i.e., based on ecological in- 

200 teractions) from genotypic evolution (i.e., based on mating-genetic interactions) speeds 

201 up diversification and approaches to the neutral scenario. Evolutionary graphs have 

202 many fascinating extensions. For example, does frequency dependence selection at ge- 

203 netic level trigger higher speed of speciation and diversity than the neutral scenario? 

204 how do gene regulatory and mating graphs interact to jump from micro to macroevo- 

205 lution? 

206 How does sexual reproduction affect evolution on graphs? Here we show that con- 

207 straining fitness according to the total number of potential matings or trophic inter- 

208 actions per individual (i.e., the genotype or phenotypic fitness model, respectively), 

209 which implies most connectivity clustered in few individuals, are a potent selection 

210 amplifier [46], and suppresses speciation rate, genetic-species diversity and species rich- 

211 ness for all the range of mutation rates and the minimum similarity values explored. 

212 This cost to diversification by common parentage factor scaling up from individuals to 

213 genetic and ecological graphs adds an additional constraint to the cost of being exces- 

214 sively abundant or rare [18] and the metabolic cost [1], thus how does the evolution 
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215 of metabolic rate interact with sexual and ecological graphs to enhance or constrain 

216 diversity at multiple biological levels and spatial scales? 

217 Most models of sympatric speciation rely on (1) intraspecific competition to drive 
21B divergence and reproductive isolation without specifying the niche or neutral nature 

219 of the interactions [13,25,40], and (2) ecological dynamics that focus on the waiting 

220 time to the first speciation [5,17]. On the other hand, neutral theory in community 

221 ecology studies patterns at the community level based on implicit modes of speciation 

222 with incipient species abundance J^ > 1 [5,18,28,30,31]. Here, despite the importance 

223 of explicit space, local adaptation and explicit prezygotic/postzygotic isolating factors 

224 to determine the mode and speed of speciation [11, 18,48,49,57,61], we link the first 

225 speciation event with the speed of speciation (i.e., mutation and fission modes of sym- 

226 patric speciation), the number of coexisting species and the genetic-species diversity 

227 in a unified framework. Note that the speed of speciation for all the parameter com- 

228 binations and models explored is extremely high. On average it is 47, 115, and 258 

229 generations to speciation, for the neutral, the phenotype and the genotype scenarios, 

230 respectively (see however [24]). If we assume a linear extrapolation from Jji (Jp) = 10'^ 

231 to 10^ inds., /i = 10"^ to 10-^ and g™" = 0.9 {Qr (Qp) ~ Qr (Qp*) ~ 0.7, see eq. 1), 

232 then the number of generations to speciation approaches to 4.7 x 10'^, 11.5 x 10^, and 

233 25.8 X 10^, which are close to the observed values in more realistic sympatric speciation 

234 models (i.e., less than 2 x 10'^ [20] and 5 x 10'^ [19] generations). 

235 Studies on food webs assume species level approaches despite the intrinsic variability 

236 in individuals [7]. In the genotypic and phenotypic fitness models only a few individuals 
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237 within each population (i.e. "hubs" ) drive reproductive rate in the context of symmetric 

238 effectiveness of ecological interactions. The expected outcome by coupling fitness with 

239 competitive and trophic asymmetry at ecological level would inevitably decrease species 

240 richness, coexistence and diversity by decreasing persistence probabilities of individuals 

241 with lower fitness. This suggest that individual variability, driven by the degree of 

242 symmetry between each pair of interacting individuals and the effectiveness of each 

243 interaction, can dramatically alter the speed of speciation, genetic-species diversity, 

244 coexistence and the structure of food webs. Note that "hubs" in networks are common 

245 but their role in inhibiting or expressing speciation and diversification at different 

246 biological levels is still unknown [59]. The need of food web data at individual level 

247 is then crucial to determine how interacting graphs at genetic and ecological levels 

248 generate the patterns of diversity and coexistence of food webs. For example, do 

249 ecological interactions depend of species or individual traits? are ecological interactions 

250 governed by a few number of individuals within each population? does neutral evolution 

251 predict the complexity and the structure of food webs? 

252 Rapid accumulation of empirical results from different biological levels suggests 

253 that ecological and evolutionary theory are undergoing a change [27,30,33]. The need 

254 to test models from first principles is now widely recognized [18,35,43,47,63]. Here 

255 we present a unified neutral model of evolution with mutation, mating with random 

256 mixing of genes, genetic-ecological drift and neutral interactions as the driving forces 

257 of diversity at multiple levels in three different scenarios. Promisingly, the huge amount 

258 of data collected and meticulously cataloged at each biological level can be used to test 
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259 neutral models from first principles in a general niche-neutral continuum multilevel 

260 framework [22, 30, 35, 43, 51, 56, 64] . 

261 4 Methods: Unifying Molecular and Ecological Evo- 

262 lution 

263 We first describe the Higgs and Derrida model of neutral molecular evolution [26] with 

264 explicit mechanisms of sympatric speciation [6]. Second, we describe the Hubbell's 

265 neutral model of biodiversity [30] with implicit speciation. Third, we highlight their 

266 similarities and link those models in the context of two initial populations that give rise 

267 to two mutualistically interacting communities. Finally we show how this framework 

268 allow us to compare the speed of speciation and the genetic-species diversity between 

269 the neutral scenario and two models of evolution at genotypic and phenotypic levels [23] 

270 using genetic and ecological graphs [46], respectively. 

271 4.1 Neutral Molecular Evolution Model 

272 Our starting point is a basic stochastic model for species formation by Higgs and 

273 Derrida (1992). This model contains three nonadaptive evolutionary forces in the 

274 sense that they are not a function of the fitness properties of the individuals: 1) 

275 neutral mutation rate (/x) in diploid and hermaphroditic individuals with equal and 

276 independent changes across any locus in a infinite genome size [38]; 2) mating with 

277 neutral mixing of genes from an hermaphroditic or two nonidentical parents and 3) 

278 genetic drift, which ensures that gene frequencies will deviate slightly from generation 
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279 to generation independent of other forces [26,47] (see section Al in the Appendix). 

280 Consider one initial completely connected and randomly mating population of size 

281 J, where individuals have the same genetic sequence. The initial genetic similarity 

282 values between any pair of individuals (g*-') in the genetic similarity matrix Q = [g*-'] 

283 are equal to 1. At each time step, an individual is chosen for death and two individuals 

284 are chosen for reproduction. Individuals have the same probability to be chosen for 

285 reproduction or for death (1/J). The viability of the offspring is constrained by g™-*", 

286 defined as the minimum genetic similarity value for postzygotic reproductive isolation 

287 (RI) two individuals i and j must satisfy for the development of fertile offspring [11, 

288 18,26,50,65]. Thus, in a randomly mating population this minimum value works as a 

289 filter generating viable offspring if and only if q^^ > g™-*". The offspring of these two 

290 individuals replace the individual that died. 

291 In this model, if the mutation rate is low (/i << 1), then the mean similarity value 

292 for Q has a solution [26, 38] 



Q* 



4J/i + l 



293 where 4J/i = 6m- The mean value arises because of a balance between mutations 

294 (which decrease the average similarity value, (Q)) and the common parentage factor 

295 which is given by the probability that two individuals have a common ancestor (which 

296 increase (Q))- Similarly, the probability that two individuals do not have a common 
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297 ancestor at stationarity is given by 



^-Q* = 7r^, (2) 



298 Results from eqs. (1) and (2) are similar to the probability that one individual is 

299 homozygous or heterozygous for one single locus under the neutral molecular theory, 

300 respectively [38]. If g™*" < Q* we will have always one species with size J . Eq. (2) 

301 represents a measure of genetic diversity in the population J . 

302 Interestingly, if g™"*" > Q* , then the initial population J is greatly perturbed by the 

303 cutoff, which implies that the genetic similarity matrix (Q) can never reach its equi- 

304 librium state. As a consequence, the initial population splits and speciation happens 

305 with the species fluctuating in the system according to demographic stochasticity [26] 

306 (see section Al in the Appendix). 

307 4.2 Neutral Theory of Biodiversity 

308 The neutral theory of biodiversity considers species instead of alleles and introduces 

309 the implicit speciation rate at the individual level (z/). The standard evolutionary 

310 metacommunity model assumes that at each time step one individual is chosen to die 

311 with probability 1/ J and is replaced by the newborn. With probability l — v-, this new 

312 individual is of the same species as its parent; with probability z/, it is an entirely new 

313 species [30,54]. 



16 



In this context, the probabihty that individual i and j in population J chosen at 
random will be of the same species is 



f - ^ '^> 



314 where the biodiversity number (Oh) is equal to Ju. 

315 In this scenario ecological drift dominates. Each individual has a percapita prob- 

316 ability to speciate at each reproduction event. This point mutation model of implicit 

317 speciation is an individual level process that leads to a proportional relationship be- 

318 tween the speciation rate of each species in a community and its abundance [14,30]. 

319 In summary, these two neutral models describe a zero-sum evolving population 

320 of J individuals with overlapping generations under demographic stochasticity. The 

321 neutral molecular model starts with a completely connected graph with mutations 

322 and mating with random mixing of genes adding variation to the new individual with 

323 explicit speciation if g™*" > Q*, The biodiversity model includes the implicit speciation 

324 parameter u. In the next section we link the implicit speciation rate (z/) with explicit 

325 mechanisms such as mutation rate (/i) and the minimum similarity value for postzygotic 

326 reproductive isolation (g™*"). Both neutral models are based on one initial population 

327 that gives rise to one community. In the next sections we describe the link between 

328 neutral molecular and biodiversity theory in the context of two initial populations that 

329 will give rise to two interacting communities. 
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330 4.3 Unified Neutral Model: Molecular and Ecological Evolu- 

331 tion 

332 The dynamics of our first two-community model has stochastic birth and death as 

333 Hubbell's model, but considers mutation (/i) and the minimum similarity value (g'"*") 

334 as explicit mechanisms of speciation at the molecular level [26] and the effectiveness of 

335 each mutualistic interaction at ecological level. We consider sexual diploid populations 

336 with overlapping generations and age independent birth and death rates. The individ- 

337 ual interactions occur in a single homogeneous patch [30,52]. The two populations can 

338 be thought of as hermaphroditic plants and dioecious pollinators, which respectively 

339 are labeled R = 1, . . . , Jr and P = 1, . . . , Jp, where Jr and Jp are the total number of 

340 individual plants and pollinators, respectively. The total number of individuals is Jm 

341 = Jr + Jp, which implies that all individuals are considered to be of reproductive age 

342 in the metacommunity. 

343 The basic model has three input parameters (i.e., the mutation rate assuming fip 

344 = fip = n, the minimum genetic similarity value assuming g™" = g™" = g"^™^ and the 

345 effectiveness of each mutualistic ecological interaction, m, assuming a fully symmetric 

346 case for all the individuals interactions) and two explicit biological levels: 1) genetic, 

347 represented as mutation (/i), mating with random mixing of genes, genetic drift and the 

348 g™"*" as in the Higgs and Derrida model already described [26,38,47], and 2) ecological 

349 level represented as ecological drift as in Hubbell's model [30] in the context of equal and 

350 symmetric competitive and mutualistic {m) ability among all interacting individuals. 

351 This is the simplest neutral scenario but the framework allows extensions to more 
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352 complicated ones. 

The speed of speciation is then governed in each plant and pollinator species by the 
mutation rate (fi), g™"*" for g™*" > Q*, and the type of sexual reproduction (i.e., facul- 
tative or obligate, mediated by the mutualistic effectiveness parameter, m). The rate 
of decay of genetic similarity of the newborn j given the similarity between the parents 
{ki{j) ,k2{j)) of the new individual j and each individual i already in the population is 
given by (see appendix): 

g^-*=__(g^iO> + ^A.,0-).)^ (4) 

353 where k^lj) can be the same than k2{j) in the hermaphroditic plant species. The time 

354 evolution of the plant and the pollinator species are governed by the generalized birth 

355 and death process with the probability of speciation in the hermaphroditic plant (z/|j) 

356 and dioecious pollinator (z/p) species k represented as: 



(5) 



-'P= .r...l .. > > Pr\ (6) 



1''' - 


Z V V pi' 2 


p — 


Ar|;(Ar^-l)Z. 2-. ^ ' 






where P^l' ^ and Pp' ^ are defined as the probabilities to produce a new species 
from two randomly chosen individuals {k^,k2) in the plant or pollinator species k, 
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respectively: 






E--i^(g*-(?'^^ + ?'^0) 




P^l'^2 

^R 


= F 


Jr-1 




t t 

^P 


= F 


Jp-1 



(7) 



where qt = ^ri^r, F[x) = 1 if x = 1 and zero otherwise. H[a) is 

if fa) = 



1 if a > 

otherwise 



Note that we have two modes of speciation. Expressions above characterize a mu- 

k' k' 
tation mode of speciation. P^^' ^ either is 1 when speciation occurs and zero otherwise. 

When the offspring of two individuals is a new individual i that cannot mate with 

any individual in the community (with i 7^ j), we have an incipient species of size 1. 

However, death events may induce the formation of new species of larger sizes. When 

there is only one individual connecting two mating networks and this happens to die, 

a new species arises. This speciation process can be called a fission speciation mode. 

In summary, our unified neutral model represents the stochastic evolution of two 

initial finite populations that give rise to two interacting communities (see Fig. 2a). 

Therefore, the interactions among individuals belonging to two different communities 

trigger the development of the ecological network as a consequence of the neutral 

dynamics at molecular and ecological levels. 
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371 4.4 Evolution on Graphs: The Genotype and the Phenotype 

372 Fitness Speciation Model 

373 Do genetic or ecological mechanisms drive the speed of speciation, genetic-species 

374 diversity and coexistence? To illuminate this question we describe two alternative 

375 models of evolution on graphs with explicit individual fitness differences. Our goal is 

376 to compare them with our unified neutral model, introduced in the last section. Fitness 

377 is defined for each individual as the reproductive probability according to the genetic 

378 similarity (i.e., genotype fitness model) or ecological affinity (i.e., phenotype fitness 

379 model) with other individuals in the same population or with individuals of the other 

380 community, respectively, but at the same time we keep neutral mutations at the genetic 

381 level symmetric. Apart from the asymmetry introduced by the different reproductive 

382 probabilities at individual level, no further asymmetry is assumed. 

383 4.4.1 The Genotype— Fitness Speciation Model 

384 Let us introduce evolution on genetic graphs as follows. In a community, individuals 

385 are labeled i = 1,2,..., J^ (Jp)- Each individual i can be described as belonging to 

386 a genetic group [45]. In each genetic group, fitness of each individual i within each 

387 species k is given by the total number of individuals j satisfying g*-' > g™-*"^ i. e., the 
383 total number of individuals each individual i can mate with. Reproductive probability 

389 of individual i within each species increases with the number of links or the number of 

390 genetically related mating partners (Fig. 2b). Thus the genetic level, using the genetic 

391 similarity among individuals, determines the speed and evolution of speciation rate and 
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392 the genetic-species diversity. Each individual i of species k is chosen for reproduction 

393 with probabihty proportional to its fitness (see appendix): 



^Nk 



SnMk 



Pi,k r, /\/r y^) 



where H{q) is 

I otherwise 

394 A^fc, Sr and M^ are the abundance of species /c, the total number of extant species in 

395 community Jr and the total number of potential mating interactions within the species 

396 k, respectively. This genotype fitness model has the following same ingredient than the 

397 unified neutral model: (1) individuals have the same probability 1/ Jr (1/Jp) to be 

398 chosen for death [60], and (2) individuals equally connected within their own species 

399 or between species are equivalent in fitness, i.e, the identity to a given species does not 

400 confer per se fitness advantage, and the following additions: (1) fitness differences are 

401 then considered only within each species according to the number of genetically related 

402 mating partners; (2) there is density dependence across species, thus rare species have 

403 higher probabilities of reproduction, and (3) we select the most connected parents with 

404 higher probability which implies that the offspring can inherits their connectance, thus 

405 increasing its reproductive probability. Evolution selects for well connected individuals. 

406 In the same way, individuals are chosen for death and reproduction in the second 

407 community. 
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408 4.4.2 The Phenotype— Fitness Speciation Model 

409 Let us now introduce evolution on ecological graphs as follows. In this last model, 

the phenotype class can be defined at the ecological level. Fitness of each individual i 

1 within each population k is associated with specific behavioral, morphological traits or 

2 phenotypic plasticity. Fitness is given by the total number of trophic links individual 

3 i of population k in one community has with j individuals belonging to populations 

4 of the other community (Fig. 2c). In this phenotype fitness model the evolution of 

5 the connectivity at the individual level within each species determines the properties 

6 of the evolving system. The reproductive probability of individual i increases with 

7 the number of trophic links but it is independent of its number of genetically related 

8 matings. Then, each individual i of species k is chosen for reproduction with probability 

9 proportional to its fitness (see appendix): 

420 where the sum until Jp means the total number of interactions of individual i with 

421 all the individuals of community Jp. rriij means that there is an interaction between 

422 individual i and j. Sr and Mk are the total number of extant species in community Jp, 

423 and the total number of mutualistic interactions among all the individuals of species k 

424 with all the individuals in community Jp, respectively. This phenotype fitness model 

425 has the same two ingredients to the neutral and genotype model: (1) individuals have 

426 the same probability l/Jp {)-/ Jp) to be chosen for death, and (2) individuals equally 
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427 connected within their own species or between species are equivalent in fitness, i.e, 

428 the identity to a given species does not confer per se fitness advantage. The model 

429 has the following additions: (1) fitness differences are then considered only within 

430 each species according to the number of trophic links with individuals of the second 

431 community, (2) there is density dependence in fitness across species, thus rare species 

432 have higher probabilities of reproduction, and (3) contrary to the genotype model, 

433 where the offspring inherits a number of potential matings from its parents, we assume 

434 that each offspring in this model starts with one trophic interaction. In the same way, 

435 individuals are chosen for death and reproduction in the second community. 
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579 6 Figure Legends 

580 • Figure 1. The Higgs and Derrida model describes the stochastic evolution of a 

581 finite population of constant size. Individuals occupy the vertex of a graph. We start 

582 with a completely connected graph with two initial populations each with Jr and Jp 

583 individuals (Fig. la). A link between each pair of individuals denotes reproductive 

584 compatibility (i.e., g*-' > g™*"). At each time step, an individual of each population is 

585 chosen for death (Fig lb). Two individuals are chosen for reproduction. Individuals 

586 have the same probability 1/ Jr {1/ Jp) to be chosen for reproduction or for death 

587 (Fig Ic). The offspring of these two individuals replace the dead individual. The 

588 newborns in Jp and Jp can be the consequence of sexual reproduction without a 

589 mutualistic interaction (i.e., facultative mutualism given by cu > 0) or a consequence of 

590 a mutualistic interaction with individual effectiveness m between the first two chosen 

591 individuals for reproduction in population Jp and Jp (i.e.. Fig. Id, given by 1 - oo). 

592 • Figure 2a represents an example of the unified neutral model. In this example 

593 each community has 5 isolated groups with different number of individuals. The most 

594 abundant groups are interacting frequently, while the rare groups do not interact among 

595 them. This neutral model is the special case of an evolving multilevel graph with fitness 

596 of each individual according to the abundance of each population. Figure 2b and 2c 

597 represent a simple scenario for the genotype and phenotype models, respectively. In the 

598 genotype scenario an individual plant and a pollinator are linked according to the total 

599 number of genetically related mating partners each individual has in its population. For 
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600 example, individuals represented with larger black circles in the Jr and Jp community 

601 have the highest number of mating links (3 in step 1). These individuals are selected 

602 as a parent of the offspring (not represented) and they are linked in step 2 (Fig. 2b). 

603 In the phenotype model an individual plant and pollinator are linked according to the 

604 total number of trophic interactions each individual has with individuals of the second 

605 community. For example, individuals with larger black circles in Jr and Jp have the 

606 highest number of trophic links (3 in step 1). These individuals are selected as the first 

607 parent of the offspring (not represented) and they are linked in step 2 (Fig. 2c). 

608 • Figure 3 represents the distribution of Spearman's rank coefficient values with Jr 

609 = Jp = 10^, /i = 10-3, and g™" equal to 0.9 {Qr (Qp) ~ Qr (Qp*) ~ 0.7, see eq. 

610 1). Fig. 3a,b,c represent the Jr community under the neutral (NUM), the phenotype 

611 (PF) and the genotype (GF) scenarios, respectively. As expected, mating and trophic 

612 graphs are not correlated during the evolution of the system under the neutral and the 

613 phenotype fitness scenarios (i.e., the distribution of Spearman's p values are close to a 

614 uniform). We randomly sampled 10^ transient values from 10 replicates with the above 

615 mentioned parameter values. Fig. 3d,e represent the individual mating (Fig. 3d) and 

616 trophic (Fig. 3e) degree ranked from the most (left) to the least (right) connected 

617 individual after 9 randomly selected transients in one replicate from the genotype 

618 fitness model. Mating and trophic degree are correlated in the genotype fitness case. 

619 The distribution of the Spearman's values is highly skewed with approximately 80% of 

620 p- values < 0.01. 

621 • Figure 4a represents speciation events as a function of the number of generations 
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622 for community Jr (community Jp not shown in the figure). Neutral, genotype and 

623 phenotype models are represented as circles, continuous and discontinuous lines, re- 

624 spectively. Each data point represents the average value after 100 replicates. We run 

625 each replicate for 10'^ generations with Jr = Jp = 10^, fi = 10"^, and g™*" equal to 0.9 

626 {Qr {Qp) ~ Qr* {Qp*) ~ 0.7, see eq. 1). Speciation events scale with the number of 

627 generations for all the three models (r^ = 0.99, red lines) with the scaling exponent 7 = 

628 0.97 (neutral), 1.03 (phenotype fitness) and 1.31 (genotype fitness). Fig. 4b represents 

629 the cumulative distribution of the number of generations to speciation for the neutral 

630 (circles), the genotype (continuous line) and the phenotype fitness models (discontin- 

631 uous line). The distributions were generated using the mean of the sorted from the 

632 smallest to the largest number of generations to speciation after 100 replicates using 

633 the same parameter values than for Fig. 4a. On average, neutral evolution generates 

634 five and twice as many speciation events as evolving genetic and ecological networks 

635 with fitness differences, respectively. 

636 • Figures 5a,b,c represent the evolution of the genetic {1—Qr), species diversity (5*6^ 

637 following eq. 10), and species richness for the neutral (circles), the genotype (continuous 

638 line) , and the phenotype fitness model (red line) with the number of generations. Jr 

639 = Jp = 10^, /i = 10-^ and g™" equal to 0.9 {Qr (Qp) ~ Qr* (Qp*) ~ 0.7, see eq. 

640 1 (results for Jp not shown are qualitatively similar). On average, neutral evolution 

641 generates twice as many genetic (Fig. 5a) and species (Fig. 5b) diversity as evolving 

642 genetic and ecological networks with fitness differences, respectively. On the other 

643 hand, on average it generates five and twice as many number of coexisting species as 

36 



644 evolving genetic and ecological networks with fitness differences, respectively (Fig. 5c). 

645 Results are the mean of 10*^ values (i.e., one value per generation) after 100 replicates 

646 with the above mentioned parameter values. 
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7 Figures 
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1 A Appendix 

2 We first describe the model at molecular level starting by the process of mutation, 

3 coancestry in overlapping generations and speciation. In the end of the first part we 

4 describe in more detail the effective number of alleles using the solution from the genetic 

5 similarity matrix in a randomly mating population (section A.l). Second we describe 

6 the effective number of species from the biodiversity theory (section A. 2). Third, we 

7 present a preliminary mathematical description of our unified neutral model by using a 

8 master equation approach (section A. 3). Fourth, we give a thorough description of our 

9 genotype (section A-4) and phenotype (section A-5) fitness speciation model. Fifth, 

10 we provide further information about how we sampled the transients and the steady 

11 state of our simulation models (section A. 6). We also have included in section A. 7 and 

12 A. 8 the table legend for the tables Al, A2, A3 and A4. Table 1 clarifies the acronyms 

13 used in the main ms. Tables 2-4 clarify the different levels and mechanisms that can 

14 be considered in a general multilevel unified model. 

15 A.l Neutral Molecular Evolution Model: Mutation, Coances- 

16 try in Overlapping Generations and Speciation 

17 In the section "Neutral Molecular Evolution Model" of main text of this paper we 

18 described briefiy the three main components of the Higgs and Derrida model: mutation 

19 (/i) , mating with random mixing of genes and genetic drift. Here we attempt to describe 

20 in more detail those components. Let fi be the average rate of mutation of the alleles 

21 existing in a diploid population. Mutation rates are equal for forward and backward 



22 mutations and across loci. We consider one initial completely connected population 

23 with all sexual identical diploids J individuals. 

24 Each individual i is represented by a sequence of A^ alleles each of which has two 

25 possible forms (+1/ — 1): (S*^, S'2...S']y), where S*^ is the -u*^ unit in the sequence of 

26 individual i [7]. The initial genetic similarity values between each pair of individuals 

27 (g*-') in the genetic similarity matrix Q = [q^^] are equal to 1. The genetic similarity 

28 between individual i and individual j is defined as 



N 
I , 



Y A^ ^ "' i-^'^^ 



^ N 

29 At each time step, an individual is chosen for death. A second individual and its 

30 partner are chosen for reproduction. Individuals have the same reproductive probability 

31 (1/J) to be chosen for reproduction or for death. If (f^ > g"^*"^ then the offspring of 

32 the two chosen parents replace (i.e., it occupies the empty site) the dead individual 

33 [4,7,11,17]. Otherwise the individuals do not mate (prezygotic RI) or their offspring 

34 is inviable or sterile (postzygotic RI) (i.e., we disregard the second individual and put 

35 it back). It is interesting to note here that this dynamics of speciation is derived from 

36 the underlying microevolutionary processes rather than postulated to follow a certain 

37 statistical distribution [4,8]. 

38 Which is the expected sequence of each new offspring? Each new individual (k) 

39 has two parents Gi{k) and 6*2 (A;). Each allele is inherited at random from one or 

40 other of the parent, thus ignoring linkage between neighboring alleles, but with a small 

41 probabihty of error determined by the mutation rate (/i). Thus if /x -^ 0, n /i ^ A, 



PiiSt - 


= S^'^'^) = ^{1 + e-^^) 


P2{St - 


4 


PsiSt - 


= S^-^'^) = ^{l + e-'n 


P^iSt - 


- -5f('')) = i(l-e-2^ 



then as n ^ oo, 

(1 _ ^)n ^ g-A (A-2) 

and the probabihty to have in the unit of the sequence S*^ the same allele than one 
of its parents is 

(A-3) 
(A-4) 
(A-5) 
(A-6) 

with the probabilities P = (Pi, Pa, P3, ^^4) satisfying < P^ < 1, i = 1, 2, ..., ^f* Pi = 
1. Given the value of the allele S'u of one of the parents, the expected value of that 
allele in the offspring is -E[>S'^] = e^'^^Su^ ■ Similarly, given the similarity between 
each parent of the new individual j and the individual i already in the population we 
update the similarity matrix q according to the following equation; 

43 We know the evolution of the overlap matrix in the limit A^ — * 00 in the infinite 

44 genome limit [7], because each pair of alleles contributing to g-^* comes with equal 

45 probability from one of the two possible combinations of the parents of j and individual 

46 i. By analyzing the similarity matrix Q at any given time it is possible to assign new 

47 individuals to a species. The analysis of the similarity matrix Q can be done by finding 



48 the isolated subgroups of individuals. In each time step we first check the k individuals 

49 that can mate with the newborn j (i.e., q^^ > g™*"). Second we check all k individuals 

50 that can mate with j but can not mate with the rest of i individuals in the population 

51 {q^^ < g™*" for all i). If this two criteria are satisfy, then we have an isolated group. 

52 The probability that two randomly selected individuals from the genetic similarity 

53 matrix Q have a common ancestor at time t is (see [1,3,9,13] for a detailed discussion 

54 about the similarities and differences with respect to monoecious and dioecious diploid 

55 populations with overlapping generations) 



Qt = e 



-if. 



^M'-i/i*-' 



(A- 



If /i << 1, then the genetic similarity matrix Q has a solution [7]: 



where 6m = 4:Jfi. This solution is identical to the inbreeding coefficient (F) in popu- 
lation genetics with non-overlapping generations [10, 16] meaning the probability that 
an individual will be homozygous. It is interesting to note that Q* and F have the 
same value despite that Q* is an average property of all loci in the sequence whereas 
F is defined for the infinite allele model in a single locus [7]. The probability that two 



individuals do not have a common ancestor is given by 

1 - Q* ^™ 



On 



(A-10) 



56 where 1 — Q* is a measure of heterozygosity or genetic diversity. As we commented 

57 in the main part of the ms. if g™*" > Q*, then the initial population J is greatly 

58 perturbed by the cutoff, the connectivity within the population decreases and the 

59 genetic similarity matrix (Q) can never reach its equilibrium state. Initial population 

60 splits and speciation happens with the species fluctuating in the system according to 

61 demographic stochasticity [7]. 

62 A. 2 Neutral Theory of Biodiversity: The Effective Number 

63 of Species 

64 Neutral theory of biodiversity discusses species instead of alleles. In this case, the 

65 probability that two individuals in population J chosen at random will be of the same 

66 species is 



ft 



V 



1 
J 



J 



ft-i 



(A-11) 



67 At equilibrium, assuming z/ is small, we have 



f' = ^ <^-'^> 



68 where 6h = Jv. Then, the effective number of species [Sg) in the J community is [14] 

Se = 0b + l (A-13) 

69 which represents a measure of species diversity and is identical to the Simpson's species 

70 diversity index [14,15]. In this scenario ecological drift dominates and speciation is an 

71 individual level process that leads to a proportional relationship between the speciation 

72 rate of each species in a community and its abundance [8]. How can we link neutral 

73 molecular and community ecology theories using evolving multilevel networks? In this 

74 section we describe in more detail the neutral unified model. 

75 A. 3 Neutral Unified Model: A Master Equation Approach 

76 with Explicit Speciation in Two Interacting Communi- 

77 ties. 

78 Our main goal here is to describe in further detail our basic neutral model. In or- 

79 der to do it, we provide a preliminary mathematical description of the evolutionary 

80 and ecological processes controlling community dynamics. An important point is that 

81 two initial populations give rise to two interacting communities, i. e., through a com- 

82 bined effect of speciation, death, and reproduction an ecological network connecting 



83 two mutualistic communities emerges. Model dynamics is controlled by three input 

84 parameters: the mutation rate (/ir = hp = fi), the minimum genetic similarity value 

85 (g™" = gp*" = g™*") at genetic level, ultimately controlling the speciation process, and 

86 the individual mutualistic effectiveness {m) assuming a fully symmetric case for all the 

87 individuals interactions at the ecological level. 

88 Our simulations consider a zero-sum birth and death stochastic individual based 

89 model in sexual diploids populations with overlapping generations and age indepen- 

90 dent birth and death rates in the context of neutral mutations and large genome size 

91 per individual (effectively infinite gene sequences) [4,7]. Individual interactions are 

92 introduced using a single and large-homogeneous patch (or metacommunity) in which 

93 there is a complete mixing and all individuals have the same chance of potentially 

94 interacting with each other [8,12]. Our model produces two mutualistically interacting 

95 communities — the resource R or plant community and the consumer or pollinator P 

96 community — but it can be easily extended to any kind of ecological interactions and 

97 to a larger number of interacting communities. 

98 The two communities assume hermaphroditic plant and dioecious pollinator indi- 

99 viduals. They are labeled R = 1, . . . , Jr, and P = 1, . . . , Jp, where Jp and Jp are 

100 the total number of individual plants and pollinators, respectively. Site size inside the 

101 patch for plants (pollinators) is defined so that each one contains one R (P) plant 

102 (pollinator) individual. Thus Jm = Jr + Jp is the total number of effective individ- 

103 uals in the metacommunity which implies that all individuals are considered in the 

104 reproductive age. These numbers are kept constant by assuming zero-sum dynamics. 
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105 An important remark is that our simulations are run much faster under the zero- 

106 sum rule, but community dynamics is easier to describe by dropping this assumption. 

107 Therefore, in what follows, for clarity and simplicity, we assume that death and non- 

108 mutualistic reproduction can take place independently in either community. Only 

109 mutualistic interactions involve simultaneously individuals from the two communities. 
no Zero-sum models are equivalent to their non zero-sum counterparts at stationarity [2]. 

111 Although we have run all our simulations under zero-sum dynamics (see our code in 

112 section A. 9), we are quite confident that our main results are robust and do not rely 

113 on the specific implementation of the zero-sum rule. 

114 Consider two randomly mating populations of size J^ and Jp where each individual 

115 has an infinitely large genome sequence subject to random neutral mutations. The 

116 initial genetic similarity values between each pair of individuals (g^ and q'p) in the 

117 matrices Qr = [q^^] and Qp = [qp] are equal to 1 and mutation rates fip and fip are 

118 equal among all individuals. Individuals have the same probability 1/Jr {1/ Jp) to 

119 be chosen for reproduction or for death. An individual of each population is chosen 

120 for death and two individuals of each population are chosen for reproduction. The 

121 offspring from the two selected individuals for reproduction of the same population 

122 replace the dying individual. Parameter uj is defined as the probability of facultative 

123 mutualistic interaction, i. e., the probability of having sexual reproduction in plants 

124 without pollinators and sexual reproduction in pollinators without plants. To be more 

125 precise, any /c-species within the plant {R) and pollinator (P) communities will undergo 

126 the following processes: 
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1. Death of an individual in the plant R community: 

2. Death of an individual in the pollinator P community: 

pk _£^ 0P (A-15) 



3. Non-mutualistic reproduction of an individual in the R community: 

0^ -^ > R^ (A-16) 



4. Non-mutualistic reproduction of an individual in the P community: 

h^.MUl-P'p) 

0^ -^ > P^ (A-17) 



5. Arising of a new species j as a result of a non-mutualistic reproduction in the R 
community: 

0^ -^ > W (A-18) 

6. Arising of a new species j as a result of a non-mutualistic reproduction in the P 
community: 

0^ -^ > P^ (A-19) 
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127 Furthermore, in case of mutualism is obligate, reproduction and ecological inter- 

128 action are considered to be coupled events involving a pair of individuals from each 

129 community. The strength of the mutualistic interaction is controlled hj 1 —u, which is 

130 defined as an individual mutualistic effectiveness, rrij^k j^k . Any time that a mutualistic 

131 event occurs, an interacting link connecting the individuals involved in the interaction 

132 appears. It is in this way that a dynamical ecological network connecting the two com- 

133 munities emerges. Notice that this reproduction-interaction coupled event can result 

134 in four different outputs: 



0^ 



™iv|=,Arfe^fl(l-^fl)^^p(l-^p) 



pk 



(A-20) 



0^ 



™Arfc,Arfc^fi(l-^p)^P^I 



R'' 
pj 



(A-21) 



0^ 



^n^^^n^^KPrMU^-Pp) 



pk 



(A-22) 



0^ 
0^ 



^Nk^,N%M^RPRMl>P'^ 



PJ 



(A-23) 
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where 0"^ (0'^)) represent empty sites, and dj^k (dj^k) and bj^k (&jv'=) death and birth 
per capita rates of plant and polhnator species k, respectively. M^ and Mp are the 
probabilities to pick up randomly two individuals that can actually mate [{i ,j ) with 
g« j > gmmj g^jYiong all available pairs in hermaphroditic plant species N^^ and pollinator 
species Np, respectively: 



2 "^ "^ 



,/ ./ 



M'p = P{^,J'\q'^ > g™"} = ^.(^. _ ^^ E E ^(^^'' -^™")' (A-25) 

i =1 j =i +1 

where H{a) 

1 if a > 
i7(a) = < 

U otherwise 



and v\ and v^ are the probabilities of mutation-induced speciation for species k in 
the plant and pollinator populations, respectively: 

4 = ^ — ^ y y py\ (A-26) 

^,fc = ^ ^y y p^i^''2 (A-27) 

tC-i — J. rCiy — /C-i ~T~ J- 
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137 where P^' ^ and Pp"' ^ are defined as the probabihties of producing an individual 

138 offspring belonging to a new species, after two parent individuals (k^jk^) have been 

139 randomly chosen from plant and pollinator k species, respectively. This new individual 

140 belongs to a new species provided it is sexually incompatible with any other individual 

141 in the respective community: 



P, 



^11^2 



R 






Jr, 



■jKjjKj 



Jp-1 



(A-28) 



(A-29) 



where qt 



-prjr, F[x] = 1 if a; = 1 and zero otherwise H(a) is 

otherwise 



143 Expressions above characterize speciation. P^^' ^ will be either 1 when speciation 

144 occurs and zero otherwise. When the offspring is a new individual i that cannot mate 

145 with any individual in the community (with i 7^ j), we have an incipient species of 

146 size 1. This a mutation- induced speciation event. However, note that there is also an 

147 alternative speciation mode. Death events may induce the formation of new species 

148 of larger sizes. When there is only one individual connecting two mating subnetworks 
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149 within the same species and this "connector" individual happens to die, the ancestral 

150 species splits into two daughter species, which implies that an additional new species 

151 arises. This is a fission- induced speciation event. In order to fully characterize this sec- 

152 ond speciation mode, we would require to count the number of "connector" individuals 

153 per species and the distribution of subnetwork sizes those individuals are connecting 

154 at any point in time. 

155 In order to write down an equation describing the temporal dynamics of the two 

156 interacting communities, we need to characterize the state of the system and all possible 

157 transitions between states at a given time. Transition rates naturally follow from the 

158 set of events we have considered above. The state of the system is defined by the 

159 community abundance vectors for the plant Nn = [N(^, N^, N^, ...Ng ] and pollinator 

160 Np = [Nf,N2,N^, ...Ng^] community, and by the two genetic similarity matrices of 

161 dimensions Jp x Jp and Jp x Jp corresponding to the two communities which control 

162 the evolving mating networks. These matrices ultimately control reproduction and 

163 speciation rates through the probabilities M^ and P^, where * stands for R and P. 

164 Notice that in the limit of large genome sizes there is no stochasticity in the similarity 

165 matrices. They are updated following the rule given by Eq. (]A-7|) after any single 

166 reproduction, extinction or speciation event. 

167 We study a fully symmetric case (m^fe ^k = rrij^k ^^ ) between plant and pollinator 

Ft^ P P '' R 

168 individuals, within species, and across species within the plant and pollinator commu- 

169 nity. In fact, for simplicity, we also assume m. .i = m,,j, ,,1.' = m = 1, i. e., the 

170 obligate mutualistic scenario, and percapita birth and death rates are assumed to take 
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the same value across species and within and across the two communities ((i|. = dp 
and 6^ = bp). This will further simplify transition rates by scaling time according to 
the birth-death temporal rate (a single time step is the time required for one death 
and one birth per community to occur). 



We define el 



[ef, . . . ,egj, where * stands for R and P and e^ is 1 ii k = i 



and otherwise. The events considered in flA-14p - flA-23p allow to define the following 
transition rates. Note these three remarks in order. First, the system can only loose 
one individual from a given species either through a death in the plant or pollinator 
community and this transition probability rate should increase linearly with the abun- 
dance of that species. Second, the encounter of two individuals for reproduction is 
a quadratic process that involve the abundance of that species squared. Third, the 
speciation process increases the dimension of the abundance community vector with 
the addition of a new component corresponding to the new incipient species. With this 
in mind, we can readily write: 



1. Death: 



T 



T 



Nn-e^Np \ Np,Np 



Np.Np-el I Np.Np 



Nh 



jf' ^^'^°^ 



j^] (A-31) 



2. Non-mutualistic reproduction: 



T 



NR + e^,Np I Nr,Np 



ujM 



{N' 



k\2 



R 



JrJr 



(1-P^) (A-32) 
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NR,Np + el I Nr,Np 



k\2 



uoM^^^^ {1 - P^) (A-33) 

JpJr 



pjp 



3. Speciation associated to non-mutualistic reproduction: 



-R 



NR + e^^^^,Np I Nr,Np 



^M^V^(^S (A-34) 



iVfl, ATp + ef^+i I Nr,Np 



ujM^^^Q^{PJ^) (A-35) 
Jpjp 



4. Mutualistic reproduction: 



T 



. (AT'^^^ 



Ar^ + ef,Arp + ef, | Nr,Np = (l-c.)M^y^ (l - P^) M 



R-'iJ 



J pj P 



I -PI 



p 
(A-36) 



5. Speciation associated to mutualistic reproduction furnishes the three remaining 
transitions (either or both species involve in the interaction undergo speciation): 



Tk \2 



NR + e^^_,„Np + e^, \ Nr,Np = (1 -a;)M^L-^ (p|) M|i ^ ' '^ ^' 



JrJi 



R 



JpJ, 



pJp 



I -PI 



(A-37) 



T 



-.R 



A^p + e« A^P + e^^+i | A^p, A^p 



Tk\2 



^R^^{i-p'k)Mi^SJ^{p 



'1 — U})M^ ^ ''"' ^^^ i\/rk' K-^" ) Ink 



JrJi 



R'JR 



J pJ p 



(A-38) 



-.R 



NR + e^^^„Np + e^^^, \ Nr,Np 



:i - uj)M: 



kiN')'^,,,,,'{N 



k \2 



R 



JrJr 



{Pr) M\ 



J pJ P 



P'. 



(A-39) 
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185 As we have done in all our simulations, if we consider only the obligate mutualistic 

1B6 scenario, cj = 0, or mutualistic effectiveness, rrij^k j^k = rrij^k j^k is equal to 1, stochas- 
1B7 tic transition rates simplify and system is described by death rates flA-30|l - flA-3ip . 



IBS mutualistic pure reproduction (lA-36p . and the speciation rates associated to obligate 
1B9 mutualistic reproduction, (lA-37P - (lA-39p . 



190 It is important to remark that the stochastic events we have considered are neglect- 

191 ing the fission speciation mode. However, all results presented in the main ms. are 

192 based on a zero-sum code which does take into account both mechanisms of speciation. 

193 Therefore, our preliminary mathematical description above provides only an approxi- 

194 mation to the actual dynamics of the system. Within this limitation, we have provided 

195 a set of transition probabilities that allow to build exact stochastic simulations [5], 

196 and a master equation that provides a general description of the time evolution of two 

197 interacting communities with explicit speciation and mutualistic interactions. These 
19B tools promise to expand our ability for quantitative analysis. However, more work is 

199 needed to fully characterize the contribution to the formation of new species through 

200 the fission-mode speciation mechanism. 

201 A. 4 Genotype Fitness Speciation Model 

202 Do genetic or ecological mechanisms drive the speed of speciation, genetic-species 

203 diversity and coexistence? We here describe in detail the individual fitness according to 

204 its reproductive probabilities. Each individual i of species k is chosen for reproduction 

205 with probability proportional to its fitness. 
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P^,k = MFi^k (A-40) 

where individual fitness is defined as: 



Thus, we write: 

\~^^k TJ ( „ij „min\ 

where A/" is a normahzation factor, N^ is the abundance of species k, and M^ is the 
total number of potential mating interactions within the species /c, which, in turn, can 
be written as: 

Nk Nk 

and H{a) 

I otherwise 
We now calculate the normalization factor by using the normalization requirement, 

i.e., by summing Pi^ across all individuals and species 1 must be obtained: 

Sr Affe 



and, then: 






Therefore, the probability of birth for each i individual is: 

Fih 



^^.^ = ^Sr ^k p. (A-46) 
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After the simplification we have: 



E^k TJ ( ^i'i „min\ 

nu= ^-^""^^ ' ^ (A-47) 

where H{a) is 

I otherwise 

206 which is the eq. 11 in the main text. N^, Sfi and M^ are the abundance of species k, 

207 the total number of extant species in community Jr and the total number of potential 

208 mating interactions within the species k, respectively. 

209 This genotype fitness model has the following same ingredient than the neutral 

210 unified model: (1) individuals have the same probability 1/Jr (l/Jp) to be chosen for 

211 death, and the following additions: (1) individuals equally connected within their own 

212 species are equivalent in fitness, i.e, the identity to a given species does not confer per 

213 se fitness advantage. Fitness differences are then considered only within each species; 

214 (2) there is a density dependence in fitness across species, thus rare species have rel- 

215 atively higher probabilities of reproduction in comparison to our unified basic neutral 

216 model described above. Note also that in this model may happen that the offspring of 

217 highly connected parents can inherit their connectance, thus increasing its reproduc- 
21B five probability. In the same way, individuals are chosen for death and reproduction 

219 in the second community. We have explored this model with the probability for death 

220 inversely proportional to the reproductive probability and the results remain qualita- 

221 tively similar to the results with the same probability to be chosen for death. In the 

222 same way, individuals are chosen for death and reproduction in the second community 
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223 Jp, thus the description of the model is equivalent. 

224 A. 5 Phenotype Fitness Speciation Model 

225 Similar to the neutral and the genotype fitness model individuals have the same prob- 

226 ability 1/ Jr {1/ Jp) to be chosen for death. Each individual i of species k is chosen for 

227 reproduction with probability proportional to their fitness, thus 






-,B _ ^j=l ""^3 



^1 = ^^^^^^ (A-48) 



where the sum until Jp means the total number of interactions of individual i with all 
individuals of community Jp. M^ is the total number of mutualistic interactions of all 
the individuals of species k with all the individuals in community Jp: 

Nk Jp 

Mfc = 5^^m,,- (A-49) 

We now use the normalization factor across all species: 

Sr. Affe 

^« = EE^5 (A-50) 

and the probability of individual i of having a newborn is: 

p.. = g (A-51) 

Finally, after the simplification we have 



^-%? <--> 
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228 which is the eq. 12 in the main text. The sum until Jp means the total number of 

229 interactions of individual i with all the individuals of community Jp. rriij means that 

230 there is an interaction between individual i and j. Sr and Mk are the total number 

231 of extant species in community J/?, and the total number of mutualistic interactions 

232 among all the individuals of species k with all the individuals in community Jp, re- 

233 spectively. This phenotype fitness model has the same two ingredients to the genotype 

234 model but working at ecological level: (1) individuals equally connected within their 

235 own species are equivalent in fitness, i.e, the identity to a given species does not con- 

236 fer per se fitness advantage. Fitness differences are then considered only within each 

237 species, (2) there is a density dependence in fitness across species, thus rare species 

238 have higher probabilities of reproduction, and (3) contrary to the genotype model, 

239 where the offspring inherits a number of potential matings from its parents, we assume 

240 that each offspring in this model starts with one trophic interaction. In the same way, 

241 individuals are chosen for death and reproduction in the second community. Similar to 

242 the genotype fitness model we have explored this model with the probability for death 

243 inversely proportional to the reproductive probability. Results remain qualitatively 

244 similar to the results with the same probability to be chosen for death. 

245 A. 6 Sampling transients and the Steady State 

246 Recent work has emphasized the importance of transient dynamics rather than long- 

247 term behavior in ecological systems [6]. In the present study we sampled transients and 

248 the steady state for the number of generations to speciation, genetic-species diversity 
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249 and species richness. This will allow us to determine if transients and the long-term 

250 behavior are similar under the neutral, the phenotype and the genotype fitness models. 

251 We thus consider both aspects of the dynamics. 

252 In all our simulations, we have assumed a fully symmetric case: rrij^k j^k = m^k j^k 

253 between plant and pollinator species, and across species within the plant and pollinator 

254 community, fn ^. ^i = m ^. ^i = m, and percapita birth and death rates have been 

255 assumed to be taken the same value across species. Furthermore, in all the replicates 

256 we have simulated the condition g™*" > Qr* (Qp*). Given Jr and Jp individuals in 

257 the initial population, a generation is an update of Jp and Jp time steps. 

258 We have explored a set of initial parameter values. Mutation rates (/i from 0.001 

259 to 0.0001), a minimum genetic similarity value to the development of viable and fertile 

260 offspring, g™*"^ from 0.75 to 0.95 in the context of a mutualistic effectiveness m = 1, and 

261 obligate mutualism, uj = 0. For the specific parameter combination of Jr = Jp = 10'^, /x 

262 = 10^^, g™*" = 0.9, m = 1, and cj = we run 100 replicates with 10^ generations each. 

263 The equilibrium value for each replicate for each community was closed to Qp (Qp) 

264 ~ Qp* (Qp*) ~ 0.7 for the neutral model. Results for all the parameter combinations 

265 explored were qualitatively similar. 
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A. 7 Table Legend of Appendix 



301 • Table Al shows the variables, parameters and acronyms used in the order that they 

302 appear in the main part of the ms. 

303 • Table A2 shows the different levels and mechanisms that can be considered in 

304 a general "Unified Neutral Model". Two general mechanisms and three levels give 

305 8 possible combinations. Note that we consider here genetic and mating as different 

306 levels. Genetic level used in the present study assumes that mutation rates are equal 

307 for forward and backward mutations and across loci. This means equal fitness among 

308 all individuals within each population at that level (i.e., /j ~ fj). Mating behavior 

309 is constrained by the minimum genetic similarity value for viable and fertile offspring 

310 (g*-' > g™*") and can be neutral as in the neutral scenario (i.e., /j ^ fj, thus all 

311 individuals within each species are equivalent) or driven by the number of genetically 

312 related matings of each individual (i.e., fi ^ fj, with explicit differences within each 

313 species) . Neutrality at ecological level assumes equivalence and symmetry in the feeding 

314 behavior {rriij = rriji = m). This neutral feeding behavior assumes competitive and 

315 mutualistic symmetric interactions of all individuals and the same percapita effect of 

316 each pollinator on each plant and viceversa. We explore here three scenarios. The first 

317 scenario represents the unified neutral model where individuals have the same fitness 

318 across all levels. This implies that each individual has the same probability to death or 

319 have descendants during the evolution of the system. This scenario is represented in the 

320 three "Neutral" conditions in the Table A2. In the phenotype and genotype models each 
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321 individual has a fitness value given by the trophic degree or number of potential matings 

322 with individuals in the second community or within the same population, respectively. 

323 Fitness is defined as the sum of the total number of j individuals each individual i 

324 interact or can mate with in each time step. The phenotype model explores the same 

325 genetic conditions than the neutral case but with different trophic interactions among 

326 individuals within each population (i.e, di ^ dj, thus the X in the continuum, Table 

327 A3). This generates differences in the reproductive probabilities (i.e., fi ^ fj). Finally, 

328 the genotype fitness model explores different mating conditions than the neutral and 

329 the phenotype models (i.e., fi ^ fj which implies the evolution of di ^ dj, thus the X 

330 at both levels. Table A4). Note that the phenotype and the genotype models end up 

331 with the same conditions but the work in opposite directions. Finally, we keep the same 

332 percapita effectiveness at ecological level in the three scenarios (i.e., rriij = rriji = m). 
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A. 8 Tables 





Table Al 


Parameter 


Meaning 


Ob 


Biodiversity number, species level 


J 


Initial population size or community size 


fi 


Mutation rate 


„min 


Genetic similarity constraint to have viable and fertile offspring 


NTME 


Neutral Theory of Molecular Evolution 


NTB 


Neutral Theory of Biodiversity 


V 


Speciation rate 


NUM 


Neutral Unified Model 


GF 


Genotypic-fitness speciation model 


PF 


Phenotypic-fitness speciation model 


r 


Incipient species size 


q^J 


Genetic similarity between individual i and j 


Q = [q''] 


Genetic similarity matrix 


RI 


Reproductive isolation 


^m 


Diversity molecular level 


Se 


Effective number of species 


Jr 


Initial population size or community size of resource/plant species 


Jp 


Initial population size or community size of pollinator species 


Jm 


Total number of individuals 


Qr = m 


Genetic similarity matrix for the plant community 


Qp = [q^] 


Genetic similarity matrix for the pollinator community 
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Level/Mechanism 


Neutral 


Niche 


Genetic {fi^i ~7^ A*) 


/. ^ /. 


/. 7^ /, 


Mating (q^^ > g™™) 


/. ^ /. 


/. 7^ fj 


Feeding 


(ij ^ dj,mij = rriji = m 


di ^ dj,mij ^ rriji 



Table A2 



Level/Mechanism 


Neutral 


Niche 


Continuum 


Genetic (/i^, ^^ /i) 


f^ ^ /. 






Mating (g^^' > g™") 




f^ ^ /. 


X 


Feeding 


rriij = rriji = "m 


di 7^ dj 


X 


Table A3: Phenotype Fitness Speciation Model 


Level/Mechanism 


Neutral 


Niche 


Continuum 


Genetic {figi^ ^^ /i) 


/. ^ fj 






Mating {q'^ > g™^") 




f^ ^ I) 


X 


Feeding 


rriij = rriji = m 


di 7^ dj 


X 



Table A4: Genotype Fitness Speciation Model 
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